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Abstract 



We study the phase diagram and critical behavior of an interacting one 
dimensional two species monomer-monomer catalytic surface reaction model 
with a reactive phase as well as two equivalent adsorbing phase where one of 
the species saturates the system. A mean field analysis including correlations 
up to triplets of sites fails to reproduce the phase diagram found by Monte 
Carlo simulations. The three phases coexist at a bicritical point whose crit- 
ical behavior is described by the even branching annihilating random walk 
universality class. This work confirms the hypothesis that the conservation 
modulo 2 of the domain walls under the dynamics at the bicritical point is the 
essential feature in producing critical behavior different from directed perco- 
lation. The interfacial fluctuations show the same universal behavior seen at 
the bicritical point in a three-species model, supporting the conjecture that 
these fluctuations are a new universal characteristic of the model. 
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I. INTRODUCTION 



Nonequilibrium statistical models with many degrees of freedom whose dynamics violate 
detailed balance arise in many areas such as biological populations, chemical reactions, fluid 
turbulence, traffic flow, and growth/deposition processes. The macroscopic behavior of these 
models can be much richer than that of systems in thermal equilibrium, showing organized 
macroscopic spatial and temporal structures like pulses or waves, and even spatiotemporal 
chaos. Even the steady state behavior of a homogeneous system without these structures 
can be far more complicated, involving for example a scale invariant steady state without 
tuning the system to a specific point. However, like their equilibrium cousins, systems at 
continuous transitions between nonequilibrium steady states show universal behavior that 
is insensitive to microscopic details and depends only on properties such as symmetries and 
conservation laws. 

One class of models that have received extensive study are those with absorbing phase 
transitions where the system changes from an active state with statistical fluctuations about 
the mean behavior to a noiseless inert state consisting of a single microscopic configuration. 
The term absorbing refers to the fact that the system cannot leave this state once it reaches 
it. Examples include directed percolation (DP) the contact process ||, auto-catalytic 
reaction models and branching annihilating random walks with odd numbers of offspring 
|5||| . Both renormalization group calculations |],|7|] and Monte Carlo simulations PHllI 
show that these models form a single universality class for a purely nonequilibrium model 
with no internal symmetry in the order parameter. 

Recently, a number of models with continuous adsorbing transitions in a universality class 
distinct from directed percolation have been studied. These models include probabilistic 
cellular automata models studied by Grassberger and coworkers [0], certain kinetic Ising 



models ||10|| , the interacting monomer-dimer model |fn}| , a three species monomer model 



with Potts-like symmetry and branching annihilating random walks with an even 

number of offspring (BAWe) p,|H|. All of these models except for the BAWe have two 
equivalent absorbing states indicating the importance of symmetry of the adsorbing state 
to the universality class. However, the universal behavior of this new class is apparently 
controlled by a dynamical conservation law. If the important dynamical variables in this class 
are defects represented by the walkers in the BAWe model and the walls between different 
saturated domains in the other models, the models have a "defect parity" conservation law 
PI where the number of defects is conserved modulo 2. Recent field theoretic work confirms 



this viewpoint [|T5[|. 

Recently, two of us [0,0 investigated a three-species monomer model with annihilation 
reactions between dissimilar species. The transitions from the reactive state to a single 
absorbing state where one species saturates the system fell in the DP universality class. 
These phase boundaries meet at bicritical points where two different absorbing states 
coexist. Because the domain walls between different domains of the two phases spawn and 
annihilate in pairs, the critical behavior fell into the BAWe universality class. We also showed 
that at the bicritical point the characteristic fluctuations of the domain walls between the 
equivalent absorbing phases was given by a new exponent. 

The present paper studies the connection between the behavior at a bicritical point 
and the presence of the BAWe critical behavior and tests the universality of the interfacial 
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fluctuations. We investigate a model introduced by Zhuo, Redner, and Park [[17] that also 
has a bicritical point at the junction of two absorbing phase transitions. We find that the 
bicritical exponents in this model also fall into the BAWe class and the interfacial exponents 
are the same as in the three species model, lending credence to the universal nature of the 
interfacial fluctuations. 

In the next section we introduce the model and present results for the mean field theory 
in two levels of approximation. Section 3 discusses our static and dynamic Monte Carlo 
simulations that yield a different phase diagram from the mean field results and produce 
critical exponents for the dynamic behavior that fall into the BAWe class. The final section 
contains our conclusions. 



II. THE MODEL 

The model we study was first introduced by Zhuo, Park and Redner [O. Two monomers, 



called A and B, adsorb at the vacant sites of a one-dimensional lattice with probabilities 
p and q, respectively, where p + q = 1. The adsorption of a monomer at a vacant site is 
affected by monomers present on neighboring sites. If either neighboring site is occupied by 
the same species as that trying to adsorb, the adsorption probability is reduced by a factor 
r < 1, mimicking the effect of a nearest-neighbor repulsive interaction. Unlike monomers 
on adjacent sites react immediately and leave the lattice, leading to a process limited only 
by adsorption. We have performed static Monte Carlo simulations that produce a phase 
diagram similar to the one found by Zhuo et al. [|17| . The diagram, displayed in Fig. [I] with 
p plotted vs. r, shows a reactive steady state (R) bordered by two equivalent saturated phases 
(labeled A and B). The transitions from the reactive phase to either of the saturated phases 
are continuous, while the transition between the saturated phases is first-order discontinuous. 



The two saturated phases meet the reactive phase at a bicritical point [16j at a critical value 
of r = r c . In the case of r = 1, the reactive region no longer exists []18| , |T9| and the only 
transition line is the first-order discontinuous line between the saturated phases. Another 
model in which the adsorption repulsion is not symmetric (only A's feel the repulsion) has 



been studied fl7j, and its critical behavior was found to be in the DP universality class. 



However, no effort was made to determine either the location of the bicritical point or the 
bicritical behavior of the model. 



III. MEAN-FIELD THEORY 

To analyze the kinetics of this model, it is useful to perform a mean-field analysis. While 
such analysis neglects long-range correlations and thus cannot be expected to properly pre- 
dict critical properties, it should properly predict the qualitative structure of the phase 
diagram, including the existence of continuous transitions and multi-critical points. The 
mean-field analysis also provides a starting point for studying the importance of such fluctu- 
ations, which become particularly important near continuous phase transitions. The mean- 
field approach we use [^0]| studies the time evolution of clusters of sites, the approximation 



coming in truncating the probabilities of observing clusters of larger size into probabilities 
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for smaller size clusters. The analysis presented below of this ID model includes clusters 
consisting of up to triplets of adjacent sites. 

At a particular time, a lattice with N sites will have Ny vacancies, the remaining sites 
being filled with Na A monomers and Nb B monomers. The density of species i = A, B, 
V is Xi = (Ni/N). We have the obvious constraint 

Xy + Xa + %B = 1, (1) 

so we can take xa and xb as independent. 

We next consider clusters of pairs of sites, where we define Nij as the number of pairs 
with species i on one site and j on the site to its right, and the pair density Xij = (Nij/N). 
Because of the immediate reaction of AB pairs we have Nab — Nba — 0. Using relations 
between the pair and site occupancies 

xa = x A a + x A v = x VA + x A a 
xb = x B b + x B v = x VB + x B b 

Xy = xav + xbv + Xw, (2) 

only two of the pair densities are independent of the x^, which we choose to be xaa and xbb- 
For the clusters of triples we define as the number of triples with species j in the 
central site, species i to the left, and species k to the right, with the corresponding densities 
Xijk — {Nijk/N). The prohibition of adjacent AB pairs and the relations between the 
numbers of clusters of triples and pairs of sites 



xaa 


= X A AV + X AAA 


= Xy AA + Xyyy 


xav 


= xaw + xav a + xavb 


= XyAV + XaAV 


xva 


= Xyy a + XAV A + XbvA 


= XyAV + XyAA 


xbb 


= xbbv + xbbb 


= xvbb + xbbb 


xbv 


= xbvv + xbv a + xbvb 


= xysv + xbbv 


xvb 


— xyyB + xavb + Xbvb 


— Xysv + XyBB 


Xyy 


— Xyyy + Xyy a + XyyB 


— Xyyy + X AW + X BW 



(3) 



gives a total of 6 independent triple densities: xaaa, xava, xbbb, xbvb, xavb, and 
xbv a- The last two are equal in a homogeneous steady state. 

We must calculate the rate that each of these densities change due to four allowed 
kinetic processes: deposition of an A without reaction, removal of an A by reaction with a 
B adsorbing next to it, deposition of a B without reaction, and removal of a B by reaction 
by an adsorbing A. These four processes are detailed in Table |. Because of the symmetry 
of the reaction rules, we will usually only present explicit formulae for half of the equations; 
the remainder can be found by interchanging A and B everywhere in the equations. 

The exact equations for the single site densities are 

P Xyyy +pr {x VVA + X AVV + Xy A v) 

- q {xwa + xaw + xvav) 

- qr(x A vB + xbv a) 



CIXa 



4 



dx B , v 

= qx vvv + qr(x VVB + x B vv + xvbv) 

— p {xvvb + %bvv + ^vbv) 

- pr (x B va + xavb) (4) 

which clearly shows the A <-> 5 symmetry noted above. 

By considering each of the possible reactions in Table Q, the equation for the pairs is 

dXAA N 

— — = pr (xyyA + ^avv + 2^^) 
at 

— g (xyyAA + %AAVV + T^^VA + T^^Aa) 

-qr {xaavb + x B vaa) (5) 

where we have denoted densities of quartets of sites in an obvious extension of our notation. 
Finally the triples obey the equations 



dx 



AAA 



dt 



dx 



AVA 



dt 



dx 



AVB 



dt 



pr (x VV AA + xaavv + xava + xavaa + xaava) 
t 1 1 ^ 

— q [XyvAAA + XAAAVV + -XAAAVA + -XavAAA) 

— qr (xaaavb + x B vaaa) 
P (xvvva + xavvv) - prXAVA 

— q (xwava + xavavv + xavava + xava) 

— qr [xavavb + x B vava) 

PXvvvb + q xavvv 

— q {xvvavb + -^xavavb) — V {xavbvv + -^xavbvb) 

-qr (xavb + x B vavb) - pr (x A vb + x A vbva) (6) 

In general, for a cluster of M sites, the adsorption processes link clusters of M sites to 
clusters of M+l sites and the reaction processes link to clusters of M+2 sites. 

To close the equations we perform a truncation. In the site approximation, we ignore any 
spatial correlations and solve just Eqn. (f|) — and its partner for i B - using Eqn. (P and 
the approximation x^k ~ XiXjXk- We find the boundary between the saturated A phase and 
the reactive phase by finding the combination of parameters p and r where the saturated 
phase fixed point loses its stability. This occurs when p = 1/(1 +r). The -B-saturated phase 
boundary is given by q = 1 — p = 1/(1 + r). These are shown as dashed lines in Fig. |2|. 
The absorbing state transitions are continuous but there is no coexistence line between the 
saturated phases, just as we found in Ref. Oj. 

In the pair approximation, we replace correlations of triples and larger clusters of sites 
by approximating them as Xijk = Xij(xjk/xj), and Xijki = Xij{xjk/xj){%ki/ x k)- Each term 
in parentheses like {xjk/xj) has the simple interpretation as the conditional probability of 
finding species k given that the site to the left is occupied by species j. We again find the 
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phase diagram by determining where the stability of the poisoned phase fixed point of the 
system of Eqns. (f|) and @ vanishes. Unlike the site approximation, solving the equations to 
determine the phase boundary was determined numerically. This procedure required some 
care because, although terms like xva and xy vanish on the phase boundary, their ratio, 
which appears in Eqns. (^) and (^) in our approximation, is nonzero. We note that, just as 
for the three species model ||13|| , Fig. ^| shows that this approximation fails to produce the 
coexistence curve between the saturated phases. This indicates that the large domains of A 
and B that appear close to the absorbing phase transitions are not accurately represented in 
the site and pair approximations. In both of these approximations, the width of the reactive 
window is proportional to (1 — r) as r approaches one. 

The triple approximation is similar in spirit to the pair approximation with higher order 
correlations being decomposed as conditional probabilities via Xijki ~ Xijk(xjki/xjk) and 
Xijkim ~ Xijk(xjki/xjk)(xkim/xki)- The numerically determined phase diagram is shown in 
Fig. and shows little difference from the pair approximation except very close to r = 1. 
Surprisingly, the triple approximation fails here to move the bicritical point away from r = 1, 
which did happen in our earlier work [jl3| on the three-species model. Very close to r = 1, 
the width of the reactive region varies as (1 — r) 3//2 in this approximation. 

We stop the mean field approach at the triple approximation because the increase in 
algebra is not compensated by an improvement of the phase diagram. The failure of this 
mean field approach to produce a realistic position for the bicritical point is linked to the 
fact that the probability of observing a long cluster of sites all filled with one species, which 
is approximated here by the probability of a smaller cluster raised to a power, will decay 
exponentially with the size of the cluster. However, the simulations presented in this paper 
and previous work |12.13| clearly show that at the bicritical point large domains of each 
species are present in the steady state, with the fundamental dynamical variables being the 
domain walls. 



IV. SIMULATIONS 

The behavior of this model at the critical lines has been studied and found to be in 



the DP universality class | l7fl , which we verified using 10 5 independent runs of up to 10 4 
time steps. What has not been studied is the behavior at the bicritical point. We use two 
types of "epidemic" simulations ||[14,21] to study the system's behavior at the bicritical 



point. The first (defect dynamics) starts with the lattice in a configuration close to one of 
the saturated phases, in which a lone vacancy is surrounded by members of one species. 
The second (interface dynamics) studies the time evolution of a minimum-width interface 
between the two saturated phases (a lone vacancy is bordered on one side by A's and by 
B's on the other side). In the defect dynamics, the simulations are over when the lattice 
becomes totally saturated. For the interface dynamics the end of the simulation is marked 
by the "collapse" of the interface back to its minimum width [12[ . 



Since the total number of possible adsorption sites is usually very small, we use a variable 
time algorithm to improve computational efficiency A list of possible adsorption sites is kept 
during the simulation, and one of these sites is picked at random for adsorption according to 
{p a } and, if applicable, r. Time is then incremented by l/riy(t), where ny(t) is the number 
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of vacancies in the lattice at that time. To avoid end effects, we always start with a lattice 
that is so large as to be considered of infinite extent. 

During the defect dynamics simulations, we measure P(t), the probability that the system 
will not become saturated at time t, ny{t), the number of vacancies in the lattice at time 
t, n Q (t), the number of species opposite to those of the initial saturated configuration at 
time t, and R 2 (t), the square size of the reactive region in a surviving run at time t. At a 
continuous phase transition as t — > oo these quantities obey power laws 

P(t) oc t~ s (R 2 (t)) oc t* (n v {t)) oc f> (n {t)) oc f ° (7) 

The exponent rj is not independent of the others, but gives a useful check on the cal- 
culation. The scaling law for r] can be understood by considering the set of all surviving 
runs. The width of the defect region grows as t z ' 2 . This region is filled with vacants and 
sites occupied by either of the species, so [B?{t)\ 1 ' 2 = c{n v + n A + n B )P(t), with c some 
constant. The factor P(t) accounts for the fact only the surviving runs contribute to R 2 (t). 
If we further assume that in this defect (which is very large for long times since z > 0) we 
have the same kind of configurations as we would see in a static simulation, we should have 
n = N^it) = iVj3 (t) . Thus we should see 

2n (t) = ^Ip^ 2 - n v (t) ~ Cl t^ 2 ~ s - c 2 f> (8) 

so we see that r\ = (z/2) — 5. 

For the interface dynamics simulations, we define P(t) to be the probability of the inter- 
face of avoiding collapse back to minimum width, with n(t) the vacancy concentration and 
R 2 (t) the square-size of the interface, both similar to the above. n (t) has no meaning or 
corresponding equivalent in this type of simulation, so it is not measured. As t — > oo these 
quantities obey power laws 

P(t)<xr 3 ', (n(t))<xf/, (R 2 (t)} oct z '. (9) 

Log-log plots of these measured quantities as a function of time are straight lines at a 
phase transition and show curvature away from transitions. A precise estimate for both the 
location of the bicritical point and the values of the exponents can be obtained by examining 
the local slopes of the curves on a log-log plot. For example, the effective exponent —S(t) is 
defined as 

-5(t) = {\n[P(t)/P(t/b)}/\nb}, (10) 

with similar definitions for r)(t), r) (t), and z(t). We usually choose 6 = 5 for our numerical 
work, with other values like b = 2 and b = 3 being used to check the results. At the bicritical 
point, the value of the effective exponent should extrapolate to the bicritical value in the 
long time limit (t^ 1 — > 0). Away from the bicritical point, the local slope will show strong 
upward or downward curvature as — > 0. 

To determine the value of both the critical exponents and their uncertainties, we use a 
technique different from that typically employed |6l,|22||. The data collected at the bicritical 
point is divided into ten data sets that are statistically independent. A linear regression is 
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then performed on each set, with the intercept the only pertinent quantity. The value for the 
exponent is then calculated as the simple mean of the ten intercepts, and the uncertainty 
is the standard error of the mean. All of the results quoted in this paper were obtained 
using this technique, which gives a more unbiased estimate of the uncertainty than previous 
methods which used the quality of the linear fit to the entire data set. 

Figure |^ shows the local slopes of the data for the defect dynamics simulations for 
p = 0.5 and r near the bicritical value. The system started near an A saturated phase, 
but the symmetry of the model dictates that we could equally as well have started near a 
B saturated phase. This data was calculated from 10 6 independent runs of up to 10 5 time 
steps at each r value. 

The analysis presented in Fig. |] gives a critical repulsion value of r c = 0.7435(15). 
We find values of 5 = 0.287(4), rj = 0.003(4), r] = 0.306(5), and z = 1.146(8) for the 
critical exponents. From these values we determine the bicritical behavior falls in the BAWe 
universality class, for which \TA\ the exponents are 5 = 0.285(2), 77 = 0.000(1), and z = 



1.141(2). The value of 7] does not agree very well with the scaling law 7] = (z/2) — 5 
derived above, but that is the result of the effect of the second term in Eq. (|), which acts 
as a large correction to the asymptotic time behavior. The agreement with the scaling law 
improves when the data for n is fit using the form (ffl) with the exponents z, 5 and 77 fixed 
and the coefficients c\ and c 2 are fitted. 

Figure ^ shows the results from the interface dynamics simulations, for three values of 
r near r c = 0.7435(15). From 10 7 independent runs of up to 10 5 time steps we find values 
of 5' = 0.700(10), rf = -0.404(9), and z' = 1.153(6). This type of simulation gives us 
information about the competition between the growth of the two saturated phases that is 
not gained in the defect dynamics type of simulations. First, notice that z', the exponent 
governing the size of the interface, equals z. Also, although 5' and rj' have different values 
than their defect dynamics cousins, the sum (5' + r/), which controls the time evolution of 
the number of vacancies in just the surviving runs, is the same as (5 + 77) within statistical 
error. This indicates that the critical spreading of the defect is insensitive as to whether it 
is bordered on each side by the same saturated phase or by two different, though equivalent, 
saturated phases. The independent dynamical exponent 5' measured here agrees with the 
value found in the three-species model 0,0], suggesting that its value is universal. 



V. SUMMARY 

We have investigated the behavior of an interacting one-dimensional adsorption-limited 
monomer-monomer model. The like-neighbor repulsion in the model leads to the presence 
of a large reactive region and a bicritical point - inaccurately termed a tricritical point 
in Ref. ]17|| - where the first-order discontinuous transition between the two symmetric 



saturated phases meets the two continuous transitions from either saturated phase to the 
reactive phase. We have used two types of epidemic simulations to study the behavior of 
the system at the bicritical point. Using the defect dynamics simulations we determine 
the bicritical behavior falls in the BAWe universality class. With the interface dynamics 
we measure a new universal number 5', related to the probability that a minimum- width 
interface between the two saturated phases will avoid collapse back to its minimum size. It 
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would be interesting to measure the exponent 5' for other models in which two symmetric 
adsorbing phases meet a reactive phase. In addition, by comparing the results from the two 
types of simulations we find that the time evolution of the number of vacancies in the lattice 
in the surviving runs is the same in both cases. This suggests that the critical spreading 
of a localized defect is insensitive to initial conditions, and the characteristic fluctuations in 
the domain walls behave differently than a defect in a homogeneous phase. 

The failure of mean field theory to describe the qualitative features of the phase diagram, 
even when clusters of triples of sites were included, points to the importance of large domains 
of each species in the vicinity of the bicritical point. The large domains, which are clearly 
apparent in static Monte Carlo simulations, are the necessary background for the simplest 
dynamical variables needed to describe the bicritical point, namely, the domain walls. The 
random walks of these walls, coupled with the fact that their number is conserved modulo 
2, are the crucial ingredients in the critical behavior at the bicritical point. 

This work was supported by the National Science Foundation under Grant No. DMR- 
9408634. 
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TABLES 



A adsorption Rate B adsorption Rate 

VVV VAV p VVV -> VBV q~ 

VVA->VAA pr VVA->VVV qr 

AW -> AAV pr AW -> VVV qr 

AVA — ► AAA pr AVA^AVV \qr 

AVA->VVA \qr 

VVB->VVV p VVB^VBB qr 

BVV->VVV p BW^BBV qr 

BVB^BW \pr BVB — > BBB qr 

BVB -► WB \pr 

AVB->AVV pr AVB^VVB qr 

BVA->VVA pr BVA->BVV qr 

TABLE I. Rates for all allowed kinetic processes. The adsorption is attempted on the middle 
site in each case 
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FIG. 1. Phase diagram determined by Monte Carlo simulation for the interacting 
monomer-monomer model. The continuous transitions from either of the saturated phases (A 
and B) to the reactive phase are shown as solid lines, while the transition between the two 
saturated phases (dot-dash line) is first-order. The three phases coexist at a bicritical point at 
( r cPc) = (0.7435, 0.5), which is shown as a heavy dot. 
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FIG. 2. Phase boundary between the A-saturated and reactive phases in the site approximation 
(solid line) and pair approximation (dashed line). The triple approximation is indistinguishable 
from the pair approximation on this figure. 
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FIG. 3. Effective exponents for the bicritical defect dynamics. The three curves correspond to 
three values of r near the bicritical point. From top to bottom they are r = 0.7420, 0.7435, 0.7450, 
with the middle curves corresponding to the bicritical values. 
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FIG. 4. Effective exponents for the bicritical interface dynamics. The three curves correspond 
to three values of r near the bicritical point. From top to bottom they are r = 0.7420, 0.7435, 0.7450, 
with the middle curve corresponding to the bicritical value. 
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